Results

library(tidyverse)
require("maps")
library(geosphere)
library(stringr)
library(rgdal)
library(caret)
library(lubridate)
library(maptools)
if (!require(ggmap)) { install.packages('ggmap'); require(ggmap) }
library(ggmap)
here_long <-  -122.3095
here_lat <- 47.6560
seattle = get_map(location = c(here_long, here_lat), zoom = 13, maptype = 'roadmap')

Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=47.656,-122.3095&zoom=13&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false

Clustering crimes by Time of Day

Time is a huge factor when discussing pedestrian safety in an error, or so we’re told. Common wisdom states that night time is more dangerous than day time, but is this even true? When do crimes get reported, and how does that change where the centers of crime are located. Here we look at a year’s worth of SPD data in order to gain some insight.

data <- read.csv('../maps-api-test/2016-2017-Clean.csv', header = TRUE)
data <- filter(data, !str_detect(Event.Clearance.Description, "HARBOR - DEBRIS, NAVIGATIONAL HAZARDS"))
ggmap(seattle) +
   geom_point(data = data, aes(x = Longitude, y = Latitude), colour = "red", alpha = 0.75)

First of, we want to make sure we use as much data as possible. Using reports for all years tract would be ideal, but that data set is too large to handle easily. Instead, we’d like to use just the past year’s worth, from November 1st, 2016, all the way to October 31, 2017. That gives us a full year’s worth of data to look at, and its far enough in the past from today that we can ensure most, if not all, incidences will be closed (and therefore included in the data set).

Before we go any further though, it is important we determine whether or not the time of year has any meaningful effect on the number of observations we have to work with. If several months have much higher crime rates than others, it may skew results of any analysis we do. With that in mind, let’s take a look at the distribution of crimes for each month in the last year:

#check frequency by month
by.month <- table(data$event_clearance_month)
ggplot(as.data.frame(by.month), 
       aes(x = Var1, y = Freq)) +
       geom_bar(stat = 'identity') +
      labs(title = 'Crimes Reported In U-District, Nov. 2016 - Oct. 2017') +
      xlab('Month') +
      ylab('Number of Crimes Reported')

As we can see, there isn’t much varience in the frequency of reported crimes in our area for the past year. We can use a Kruskal-Wallis to test the independence of Month and Number of Crimes, which should indicate whether there is a relationship between them or not:

k <- kruskal.test(Freq ~ Var1, by.month)
k

    Kruskal-Wallis rank sum test

data:  Freq by Var1
Kruskal-Wallis chi-squared = 11, df = 11, p-value = 0.4433

Since our null hypothesis, that the count of crimes for each month is consistent, was given a p-value of 0.4433, we can confidently reject it and state that there is no dependence between the month of year and the number of crimes reported during it.

Next, we can look at the distribution of crimes across the categories we’ve defined. If there’s no variation there, we can continue to look at the data set as a whole, but if there’s significant variation, we’ll need to handle things a bit more carefully. Here we see a histogram of the Event Clearance Descriptions (what the reported crime was classified as in the SPD’s computer system.):

freq_by_desc <- table(droplevels(data$Event.Clearance.Description))
ggplot(as.data.frame(freq_by_desc), 
       aes(x = Var1, y = Freq)) +
       geom_bar(stat = 'identity') +
      xlab('Number of Crimes') +
      ylab('Crime Description') +
    coord_flip()

We see here that there is significant differences in the types of crimes that are reported, with Crisis Complaints comprising a large number of them. This could simply be due to a large number of mental health crisises in the U-District. More likely, however, officers are unsure of what to classify a crime as in an incident report and they simply choose a catch all category that comes closest to describing what happened. Regardless, this uneven distribution means that when it comes time to perform clustering, we’ll need to be careful to account for the significant weight these crimes will inflict upon the cluster centers.

#ggmap(seattle) +
#  geom_point(data = data, aes(x = Longitude, y = Latitude, group = Event.Clearance.Description, color = #Event.Clearance.Description), alpha = 0.5, size = 10) +
#  facet_wrap(~ Event.Clearance.Description) +
#  theme(axis.ticks = element_blank(), 
#        axis.text.x = element_blank(),
#        axis.text.y = element_blank(),
#        strip.text = element_text(size=50),
#        legend.position = "none"
#        )

These are some function that will help facilitate running k means clustering

#some useful functions for performing clustering
#extract the lat and long from a dataframe, and run kmeans on it
# x = one of our dataframes
# clusters = how many centers you want kmeans to work with when clustering
fit.clusters <- function(x, clusters) {
  # selecting just ID and location data
  df_loc <- x %>% dplyr::select(CAD.CDW.ID, Latitude, Longitude)
  # fitting model
  fit <- kmeans(df_loc, clusters)
  fit$centers # look at cluster sizes and means. want clusters to be about equal size
  return(fit)
}
#make a plot that will tell you how many clusters might work for a given dataframe
# x = a dataframe
# max = the maximum number of clusters you want to try
find.num.clusters <- function(x, max) {
  if(max > nrow(x)) { stop('Cannot fit more clusters than there are rows in dataframe')}
  df_loc <- x %>% dplyr::select(CAD.CDW.ID, Latitude, Longitude)
  wss = c()
  for (i in 1:max) {
    wss[i] <- sum(kmeans(df_loc, centers=i)$withinss)
  }
  plot(1:max, wss, type="b", xlab="Number of Clusters",
     ylab="Within groups sum of squares")
}
#plot the number of observations in each cluster
# x = a fit object returned from kmeans() or the fit.clusters() function above
plot.cluster.sizes <- function(x) {
  cluster.size <- data.frame(data.frame('clusters' = 1:nrow(x$centers), x$size))
  ggplot(data = cluster.size, aes(x = clusters, y = x.size)) +
  geom_bar(stat = 'identity')
}

We can bucket the data into 6 time frames to look at how reported crime changes throughout the day:

morning <- filter(data, 6 <= at_scene_time_hr, at_scene_time_hr < 10 )
mid.day <-  filter(data, 10 <= at_scene_time_hr, at_scene_time_hr < 14 )
afternoon <-  filter(data, 14 <= at_scene_time_hr, at_scene_time_hr < 18 )
evening <-  filter(data, 18 <= at_scene_time_hr, at_scene_time_hr < 22 )
night <-  filter(data, 22 <= at_scene_time_hr | at_scene_time_hr < 2 )
early.morning <-  filter(data, 2 <= at_scene_time_hr, at_scene_time_hr < 6 )
ggmap(seattle) +
  geom_point(data = morning, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Morning')

ggmap(seattle) +
  geom_point(data = mid.day, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Mid Day')

ggmap(seattle) +
  geom_point(data = afternoon, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Afternoon')

ggmap(seattle) +
  geom_point(data = evening, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Evening')

ggmap(seattle) +
  geom_point(data = night, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Night')

ggmap(seattle) +
  geom_point(data = early.morning, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Early Morning')

We can see that crimes reported are generally normally distributed around the afternoon. This would suggest that the highest rate of crime is during the day, or that there are less people reporting crimes at night. Which one is true isn’t possible to infer from this data.

lengths <- c(nrow(morning), nrow(mid.day), nrow(afternoon), nrow(evening), nrow(night), nrow(early.morning))
names <- c('Morning\n6:00 - 9:59', 'Mid-day\n10:00 - 1:59', 'Afternoon\n2:00 - 5:59', 'Evening\n6:00 - 9:59', 'Night\n10:00 - 1:59', 'Early Morning\n2:00 - 5:59')
by.tod <- data.frame('TOD' = names, 'Count.Crimes' = lengths)
by.tod$TOD = factor(by.tod$TOD, levels = by.tod$TOD)
ggplot(by.tod, aes(x = TOD, y = Count.Crimes)) +
  geom_histogram(stat = 'identity')

Ignoring unknown parameters: binwidth, bins, pad

Clustering reveals that the average reported crime doesn’t stay too far from the Ave. Some time frames have somewhat lower densities.

fit <- fit.clusters(morning, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

plot.cluster.sizes(fit)

fit <- fit.clusters(mid.day, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

plot.cluster.sizes(fit)

fit <- fit.clusters(afternoon, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

plot.cluster.sizes(fit)

fit <- fit.clusters(evening, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

plot.cluster.sizes(fit)

fit <- fit.clusters(night, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

plot.cluster.sizes(fit)

The fact that Crisis Complaints outnumber every other description and its not very specific is skewing our analysis. Let’s remove those crimes and try again.

#take out general crisis complaint - general
data <- filter(data, Event.Clearance.Description != 'CRISIS COMPLAINT - GENERAL')
morning <- filter(data, 6 <= at_scene_time_hr, at_scene_time_hr < 10 )
mid.day <-  filter(data, 10 <= at_scene_time_hr, at_scene_time_hr < 14 )
afternoon <-  filter(data, 14 <= at_scene_time_hr, at_scene_time_hr < 18 )
evening <-  filter(data, 18 <= at_scene_time_hr, at_scene_time_hr < 22 )
night <-  filter(data, 22 <= at_scene_time_hr | at_scene_time_hr < 2 )
early.morning <-  filter(data, 2 <= at_scene_time_hr, at_scene_time_hr < 6 )
ggmap(seattle) +
  geom_point(data = morning, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Morning')

ggmap(seattle) +
  geom_point(data = mid.day, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Mid Day')

ggmap(seattle) +
  geom_point(data = afternoon, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Afternoon')

ggmap(seattle) +
  geom_point(data = evening, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Evening')

ggmap(seattle) +
  geom_point(data = night, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Night')

ggmap(seattle) +
  geom_point(data = early.morning, aes(x = Longitude, y = Latitude), alpha = 0.5) + 
  labs(title = 'Early Morning')

lengths <- c(nrow(morning), nrow(mid.day), nrow(afternoon), nrow(evening), nrow(night), nrow(early.morning))
names <- c('Morning\n6:00 - 9:59', 'Mid-day\n10:00 - 1:59', 'Afternoon\n2:00 - 5:59', 'Evening\n6:00 - 9:59', 'Night\n10:00 - 1:59', 'Early Morning\n2:00 - 5:59')
by.tod <- data.frame('TOD' = names, 'Count.Crimes' = lengths)
by.tod$TOD = factor(by.tod$TOD, levels = by.tod$TOD)
ggplot(by.tod, aes(x = TOD, y = Count.Crimes)) +
  geom_histogram(stat = 'identity')

Ignoring unknown parameters: binwidth, bins, pad

# find the mode of numeric/character data
Mode <- function(x) {
  ux <- unique(x)
  tab <- tabulate(match(x, ux)); ux[tab == max(tab)]
}
tod.mean <- mean(data$at_scene_time_hr)
tod.med <- median(data$at_scene_time_hr)
tod.mean

[1] 13.6691

tod.med

[1] 14

Mode(data$at_scene_time_hr)

[1] 17

#What is the most common crime committed at each period?
Mode(morning$Event.Clearance.Description)

[1] HAZARDS
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(mid.day$Event.Clearance.Description)

[1] HARASSMENT, THREATS
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(afternoon$Event.Clearance.Description)

[1] HARASSMENT, THREATS
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(evening$Event.Clearance.Description)

[1] HARASSMENT, THREATS
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(night$Event.Clearance.Description)

[1] HARASSMENT, THREATS
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(early.morning$Event.Clearance.Description)

[1] STRONG ARM ROBBERY
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

fit <- fit.clusters(morning, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

plot.cluster.sizes(fit)

fit <- fit.clusters(mid.day, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

plot.cluster.sizes(fit)

fit <- fit.clusters(afternoon, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

plot.cluster.sizes(fit)

fit <- fit.clusters(evening, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

plot.cluster.sizes(fit)

fit <- fit.clusters(night, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

plot.cluster.sizes(fit)